home *** CD-ROM | disk | FTP | other *** search
/ Developer CD Series 1998 April: Mac OS SDK / Dev.CD Apr 98 SDK1.toast / Development Kits (Disc 1) / QuickDraw 3D / Samples / SampleCode / Plug-in - WireFrame Renderer / SR_Math.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-08-14  |  22.5 KB  |  844 lines  |  [TEXT/MPS ]
  1. /******************************************************************************
  2.  **                                                                             **
  3.  **     Module:        SR_Math.c                                                 **
  4.  **                                                                          **
  5.  **                                                                          **
  6.  **     Purpose:     Various mathematical utility functions                      **
  7.  **                                                                          **
  8.  **                                                                          **
  9.  **                                                                          **
  10.  **     Copyright (C) 1996 Apple Computer, Inc.  All rights reserved.          **
  11.  **                                                                          **
  12.  **                                                                          **
  13.  *****************************************************************************/
  14. #include <assert.h>
  15. #include <math.h>
  16.  
  17. #include "QD3D.h"
  18. #include "QD3DMath.h"
  19. #include "QD3DGeometry.h"
  20.  
  21. #include "SR.h"
  22. #include "SR_Math.h"
  23.  
  24.  
  25. /*===========================================================================*\
  26.  *
  27.  *    Routine:    SRTriangle_GetNormal()
  28.  *
  29.  *    Comments:    Compute the face normal for a triangle. If the triangle
  30.  *                has a normal in its attribute set, use that; otherwise,
  31.  *                compute the normal from the vertices.
  32.  *
  33. \*===========================================================================*/
  34.  
  35. TQ3Status SRTriangle_GetNormal(
  36.     TQ3TriangleData        *triangleData,
  37.     TQ3Vector3D            *normal)
  38. {
  39.     assert(triangleData != NULL);
  40.         
  41.     /*
  42.      *  Check for face normal for triangle. Use it if it's there, but assume
  43.      *  it's normalized, as per the rule.
  44.      */
  45.     if ((triangleData->triangleAttributeSet != NULL) && 
  46.         (Q3AttributeSet_Get( 
  47.             triangleData->triangleAttributeSet, 
  48.             kQ3AttributeTypeNormal, 
  49.             normal) == kQ3Success)) {
  50.     } else {
  51.         /* 
  52.          *  Otherwise, compute plane normal, and normalize.
  53.          */
  54.         Q3Point3D_CrossProductTri(
  55.             &triangleData->vertices[0].point,
  56.             &triangleData->vertices[1].point,
  57.             &triangleData->vertices[2].point,
  58.             normal);
  59.         Q3Vector3D_Normalize(normal, normal);
  60.     }
  61.     
  62.     return (kQ3Success);    
  63. }
  64.  
  65.  
  66. /*===========================================================================*\
  67.  *
  68.  *    Routine:    SRMatrix_LUDecomposeSingular3x3()
  69.  *
  70.  *    Comments:    Decompose a singular 3x3 matrix "A" into PA = LU.
  71.  *
  72.  *                P is the permutation matrix for the row exchanges required
  73.  *                    by pivoting.
  74.  *                L is a lower triangular matrix with ones on the diagonal.
  75.  *                U is an upper echelon matrix.  For nonsingular matrices,
  76.  *                    it is upper triangular.
  77.  *
  78.  *                Upon return, this routine returns LU in place of the original
  79.  *                matrix.  U occupies the upper triangular elements.  L occupies
  80.  *                the elements below the diagonal.  P is represented in "short-
  81.  *                hand" form by the three-element vector indexPtr.  If the
  82.  *                number of row exchanges is even, then *rowInterchanges is 1.0,
  83.  *                otherwise, -1.0.  The rank of the matrix is returned in *rank.
  84.  *
  85.  *                To fully understand this routine, read Chapters 1-3 of 
  86.  *                Gilbert Strang's book, Linear Algebra and Its Applications.
  87.  *
  88. \*===========================================================================*/
  89.  
  90. /* 
  91.  *  Macros for supporting SRMatrix_LUDecomposeSingular3x3 
  92.  */
  93. #define    FLOAT_TOLERANCE        1.0e-5
  94. #define EQUALS_ZERO(a,e)     (FABS(a) < (e))
  95. #define SWAP(f,g)             \
  96.     {                         \
  97.         float t;             \
  98.                             \
  99.         t = (f);             \
  100.         (f) = (g);             \
  101.         (g) = t;             \
  102.     }
  103.  
  104. void SRMatrix_LUDecomposeSingular3x3(
  105.     float     matrix[3][3],
  106.     long     indexPtr[3],
  107.     float     *rowInterchanges,
  108.     long    *rank)
  109. {
  110.     long    i, j;
  111.     float    f;
  112.     float    pivot;
  113.     long    iPivot;
  114.     float    max[3];
  115.     
  116.     /* Initialize parameters */
  117.     *rowInterchanges = 1.0;
  118.     for (i = 0; i < 3; i++) {
  119.         indexPtr[i] = i;
  120.     }
  121.     *rank = 3;
  122.     
  123.     /* Determine magnitude of largest element in each column */
  124.     for (j = 0; j < 3; j++) {
  125.         max[j] = 0.0;
  126.         for (i = 0; i < 3; i++) {
  127.             if ((f = FABS(matrix[i][j])) > max[j]) {
  128.                 max[j] = f;
  129.             }
  130.         }
  131.     }
  132.     
  133.     /*
  134.      * Initial knowledge about LU decomposition:
  135.      *
  136.      *            |1 0 0|            |? ? ?|
  137.      *        L = |? 1 0|        U = |? ? ?|
  138.      *            |? ? 1|            |? ? ?|
  139.      *
  140.      *    ? = yet to be determined
  141.      *    * = determined final value
  142.      *  X = nonzero pivot
  143.      */
  144.      
  145.     /* Find the row with the largest pivot in the first column */
  146.     for (i = 0, pivot = 0.0, iPivot = 0; i < 3; i++) {
  147.         f = matrix[i][0];
  148.         if (FABS(f) > pivot) {
  149.             pivot = f;
  150.             iPivot = i;
  151.         }
  152.     }
  153.     
  154.     /* Record permutation of rows */
  155.     if (iPivot != 0) {
  156.         /* Swap row with largest pivot with first row */
  157.         for (j = 0; j < 3; j++) {
  158.             SWAP(matrix[0][j], matrix[iPivot][j]);
  159.         }
  160.         SWAP(indexPtr[0],indexPtr[iPivot]);
  161.         *rowInterchanges = -(*rowInterchanges);
  162.     }
  163.     
  164.     if ((max[0] == 0.0) || 
  165.          EQUALS_ZERO(matrix[0][0] / max[0], FLOAT_TOLERANCE)) {
  166.         /*
  167.          * First column of U is zero:
  168.          *
  169.          *            |1 0 0|            |0 ? ?|
  170.          *        L = |? 1 0|        U = |0 ? ?|
  171.          *            |? ? 1|            |0 ? ?|
  172.          */
  173.         (*rank)--;
  174.         matrix[0][0] = 0.0;
  175.         
  176.         /* Find the row with the largest pivot in the second column */
  177.         for (i = 0, pivot = 0.0, iPivot = 0; i < 3; i++) {
  178.             f = matrix[i][1];
  179.             if (FABS(f) > pivot) {
  180.                 pivot = f;
  181.                 iPivot = i;
  182.             }
  183.         }
  184.  
  185.         /* Record permutation of rows */
  186.         if (iPivot != 0) {
  187.             /* Swap row with largest pivot with first row */
  188.             for (j = 0; j < 3; j++) {
  189.                 SWAP(matrix[0][j], matrix[iPivot][j]);
  190.             }
  191.             SWAP(indexPtr[0],indexPtr[iPivot]);
  192.             *rowInterchanges = -(*rowInterchanges);
  193.         }
  194.  
  195.         if ((max[1] == 0.0) || 
  196.              EQUALS_ZERO(matrix[0][1] / max[1], FLOAT_TOLERANCE)) {
  197.             /*
  198.              * Second column of U is zero:
  199.              *
  200.              *            |1 0 0|            |0 0 ?|
  201.              *        L = |? 1 0|        U = |0 0 ?|
  202.              *            |? ? 1|            |0 0 ?|
  203.              */
  204.             (*rank)--;
  205.             matrix[0][1] = 0.0;
  206.             matrix[1][1] = 0.0;
  207.             
  208.             /* Find the row with the largest pivot in the third column */
  209.             for (i = 0, pivot = 0.0, iPivot = 0; i < 3; i++) {
  210.                 f = matrix[i][2];
  211.                 if (FABS(f) > pivot) {
  212.                     pivot = f;
  213.                     iPivot = i;
  214.                 }
  215.             }
  216.     
  217.             /* Record permutation of rows */
  218.             if (iPivot != 0) {
  219.                 /* Swap row with largest pivot with first row */
  220.                 for (j = 0; j < 3; j++) {
  221.                     SWAP(matrix[0][j], matrix[iPivot][j]);
  222.                 }
  223.                 SWAP(indexPtr[0],indexPtr[iPivot]);
  224.                 *rowInterchanges = -(*rowInterchanges);
  225.             }
  226.  
  227.             if ((max[2] == 0.0) || 
  228.                  EQUALS_ZERO(matrix[0][2] / max[2], FLOAT_TOLERANCE)) {
  229.                 /*
  230.                  * Third column of U is zero:
  231.                  *
  232.                  *            |1 0 0|            |0 0 0|
  233.                  *        L = |0 1 0|        U = |0 0 0|
  234.                  *            |0 0 1|            |0 0 0|
  235.                  */
  236.                 (*rank)--;
  237.                 matrix[0][2] = 0.0;
  238.                 matrix[1][0] = 0.0;
  239.                 matrix[1][2] = 0.0;
  240.                 matrix[2][0] = 0.0;
  241.                 matrix[2][1] = 0.0;
  242.                 matrix[2][2] = 0.0;
  243.             } else {
  244.                 /*
  245.                  * Third column of U has a nonzero pivot:
  246.                  * Need to calculate multiplier of L to
  247.                  * eliminate third column of of U:
  248.                  *
  249.                  *            |1 0 0|            |0 0 X|
  250.                  *        L = |* 1 0|        U = |0 0 0|
  251.                  *            |* 0 1|            |0 0 0|
  252.                  */
  253.                 matrix[1][0] = matrix[1][2] / matrix[0][2];
  254.                 matrix[1][2] = 0.0;
  255.                 matrix[2][0] = matrix[2][2] / matrix[0][2];
  256.                 matrix[2][1] = 0.0;
  257.                 matrix[2][2] = 0.0;
  258.             }
  259.         } else {
  260.             /*
  261.              * Second column of U has a nonzero pivot:
  262.              * Need to calculate multiplier of L to
  263.              * eliminate second column of of U:
  264.              *
  265.              *            |1 0 0|            |0 X *|
  266.              *        L = |* 1 0|        U = |0 0 ?|
  267.              *            |* ? 1|            |0 0 ?|
  268.              */
  269.             f = matrix[1][1] / matrix[0][1];
  270.             matrix[1][0] = f;
  271.             matrix[1][1] = 0.0;
  272.             matrix[1][2] = matrix[1][2] - f * matrix[0][2];
  273.              
  274.             f = matrix[2][1] / matrix[0][1];
  275.             matrix[2][0] = f;
  276.             matrix[2][1] = 0.0;
  277.             matrix[2][2] = matrix[2][2] - f * matrix[0][2];
  278.  
  279.             /* Find the row with the largest pivot in the third column */
  280.             if (FABS(matrix[1][2]) < FABS(matrix[2][2])) {
  281.                 iPivot = 2;
  282.                 /* Swap third row (possessing largest pivot) with second row */
  283.                 for (j = 0; j < 3; j++) {
  284.                     SWAP(matrix[1][j], matrix[iPivot][j]);
  285.                 }
  286.                 SWAP(indexPtr[1], indexPtr[iPivot]);
  287.                 *rowInterchanges = -(*rowInterchanges);
  288.             } else {
  289.                 iPivot = 1;
  290.             }
  291.             
  292.             if ((max[2] == 0.0) || 
  293.                  EQUALS_ZERO(matrix[1][2] / max[2], FLOAT_TOLERANCE)) {
  294.                 /*
  295.                  * Third column of U has no nonzero pivot.
  296.                  *
  297.                  *            |1 0 0|            |0 X *|
  298.                  *        L = |* 1 0|        U = |0 0 0|
  299.                  *            |* 0 1|            |0 0 0|
  300.                  */
  301.                 (*rank)--;
  302.                 matrix[1][2] = 0.0;
  303.                 matrix[2][1] = 0.0;
  304.                 matrix[2][2] = 0.0;
  305.             } else {
  306.                 /*
  307.                  * Third column of U has a nonzero pivot.
  308.                  * Need to calculate multiplier of L to
  309.                  * eliminate third column of of U:
  310.                  *
  311.                  *            |1 0 0|            |0 X *|
  312.                  *        L = |* 1 0|        U = |0 0 X|
  313.                  *            |* * 1|            |0 0 0|
  314.                  */
  315.                 matrix[2][1] = matrix[2][2] / matrix[1][2];;
  316.                 matrix[2][2] = 0.0;
  317.             }
  318.         }
  319.     } else {
  320.         /*
  321.          * First column of U has a nonzero pivot.
  322.          * Need to calculate multipliers of L to
  323.          * eliminate first column of of U:
  324.          *
  325.          *            |1 0 0|            |X * *|
  326.          *        L = |* 1 0|        U = |0 ? ?|
  327.          *            |* ? 1|            |0 ? ?|
  328.          */
  329.         f = matrix[1][0] / matrix[0][0];
  330.         matrix[1][0] = f;
  331.         matrix[1][1] = matrix[1][1] - f * matrix[0][1];
  332.         matrix[1][2] = matrix[1][2] - f * matrix[0][2];
  333.          
  334.         f = matrix[2][0] / matrix[0][0];
  335.         matrix[2][0] = f;
  336.         matrix[2][1] = matrix[2][1] - f * matrix[0][1];
  337.         matrix[2][2] = matrix[2][2] - f * matrix[0][2];
  338.  
  339.         /* Find the row with the largest pivot in the second column */
  340.         if (FABS(matrix[1][1]) < FABS(matrix[2][1])) {
  341.             iPivot = 2;
  342.             /* Swap third row (possessing largest pivot) with second row */
  343.             for (j = 0; j < 3; j++) {
  344.                 SWAP(matrix[1][j], matrix[iPivot][j]);
  345.             }
  346.             SWAP(indexPtr[1], indexPtr[iPivot]);
  347.             *rowInterchanges = -(*rowInterchanges);
  348.         } else {
  349.             iPivot = 1;
  350.         }
  351.         
  352.         if ((max[1] == 0.0) || 
  353.              EQUALS_ZERO(matrix[1][1] / max[1], FLOAT_TOLERANCE)) {
  354.             /*
  355.              * Second column of U has no nonzero pivot.
  356.              *
  357.              *            |1 0 0|            |X * *|
  358.              *        L = |* 1 0|        U = |0 0 ?|
  359.              *            |* ? 1|            |0 0 ?|
  360.              */
  361.             (*rank)--;
  362.             matrix[1][1] = 0.0;
  363.     
  364.             /* Find the row with the largest pivot in the third column */
  365.             if (FABS(matrix[1][2]) < FABS(matrix[2][2])) {
  366.                 iPivot = 2;
  367.                 /* Swap third row (possessing largest pivot) with second row */
  368.                 for (j = 0; j < 3; j++) {
  369.                     SWAP(matrix[1][j], matrix[iPivot][j]);
  370.                 }
  371.                 SWAP(indexPtr[1], indexPtr[iPivot]);
  372.                 *rowInterchanges = -(*rowInterchanges);
  373.             } else {
  374.                 iPivot = 1;
  375.             }
  376.             
  377.             if ((max[2] == 0.0) || 
  378.                  EQUALS_ZERO(matrix[1][2] / max[2], FLOAT_TOLERANCE)) {
  379.                 /*
  380.                  * Third column of U has no nonzero pivot.
  381.                  *
  382.                  *            |1 0 0|            |X * *|
  383.                  *        L = |* 1 0|        U = |0 0 0|
  384.                  *            |* 0 1|            |0 0 0|
  385.                  */
  386.                 (*rank)--;
  387.                 matrix[1][2] = 0.0;
  388.                 matrix[2][1] = 0.0;
  389.                 matrix[2][2] = 0.0;
  390.             } else {
  391.                 /*
  392.                  * Third column of U has a nonzero pivot.
  393.                  * Need to calculate multiplier of L to
  394.                  * eliminate third column of of U:
  395.                  *
  396.                  *            |1 0 0|            |X * *|
  397.                  *        L = |* 1 0|        U = |0 0 X|
  398.                  *            |* * 1|            |0 0 0|
  399.                  */
  400.                 matrix[2][1] = matrix[2][2] / matrix[1][2];;
  401.                 matrix[2][2] = 0.0;
  402.             }
  403.         } else {
  404.             /*
  405.              * Second column of U has a nonzero pivot.
  406.              * Need to calculate multipliers of L to
  407.              * eliminate second column of of U:
  408.              *
  409.              *            |1 0 0|            |X * *|
  410.              *        L = |* 1 0|        U = |0 X *|
  411.              *            |* * 1|            |0 0 ?|
  412.              */
  413.             f = matrix[2][1] / matrix[1][1];
  414.             matrix[2][1] = f;
  415.             matrix[2][2] = matrix[2][2] - f * matrix[1][2];
  416.             
  417.             if ((max[2] == 0.0) || 
  418.                  EQUALS_ZERO(matrix[2][2] / max[2], FLOAT_TOLERANCE)) {
  419.                 /*
  420.                  * Third column of U has no nonzero pivot.
  421.                  *
  422.                  *            |1 0 0|            |X * *|
  423.                  *        L = |* 1 0|        U = |0 X *|
  424.                  *            |* * 1|            |0 0 0|
  425.                  */
  426.                 (*rank)--;
  427.                 matrix[2][2] = 0.0;
  428.             } else {
  429.                 /*
  430.                  * Third column of U has a nonzero pivot.
  431.                  * The matrix has full rank.
  432.                  *
  433.                  *            |1 0 0|            |X * *|
  434.                  *        L = |* 1 0|        U = |0 X *|
  435.                  *            |* * 1|            |0 0 X|
  436.                  */
  437.             }
  438.         }
  439.     }
  440. }
  441.  
  442.  
  443. /*===========================================================================*\
  444.  *
  445.  *    Routine:    SRMatrix_ComputeFlatLand()
  446.  *
  447.  *    Comments:    Computes the eye vector in local coordinates and the plane
  448.  *                equation in world coordinates when the upper-left 3x3 submatrix
  449.  *                of the local-to-world matrix has a rank of 2.  The eye vector
  450.  *                in local coordinates can be used for culling, but should not
  451.  *                be used for lighting because the local-to-world matrix does
  452.  *                not preserve lengths.
  453.  *
  454.  *    Input:        matrix:                the local-to-world matrix
  455.  *                luDecomp:            the LU decomposition of the upper-left 
  456.  *                                    3x3 submatrix of the local-to-world matrix
  457.  *                parallelProjection:    true for parallel projections, false for
  458.  *                                    perspective projections
  459.  *                eyeVectorWC:        the eye vector in world coordinates
  460.  *                                    (required only for parallel projections)
  461.  *                eyePointWC:            the eye point in world coordinates
  462.  *                                    (required only for perspective projections)
  463.  *
  464.  *    Output:        eyeVectorLC:        the eye vector in local coordinates
  465.  *                flatWorld:            the plane equation in world coordinates
  466.  *                                    to which all geometry is transformed with
  467.  *                                    the normal in the same hemisphere as the 
  468.  *                                    eye
  469.  *
  470.  *    Returns:    kQ3Success:            computation completed successfully
  471.  *                kQ3Failure:            computation could not be completed 
  472.  *                                    successfully
  473.  *
  474.  *                To fully understand this routine, read Chapters 1-3 of 
  475.  *                Gilbert Strang's book, Linear Algebra and Its Applications.
  476.  *
  477. \*===========================================================================*/
  478.  
  479. TQ3Status SRMatrix_ComputeFlatLand(
  480.     TQ3Matrix4x4        *matrix,
  481.     TQ3Matrix3x3        *luDecomp,
  482.     TQ3Boolean            parallelProjection,
  483.     TSRVector4D            *eyeVectorWC,
  484.     TQ3RationalPoint4D    *eyePointWC,
  485.     TSRVector4D            *eyeVectorLC,
  486.     TQ3PlaneEquation    *flatWorld)
  487. {
  488.     long            pivot0 = -1;
  489.     long            pivot1 = -1;
  490.     long            j;
  491.     TQ3Vector3D        v0;
  492.     TQ3Vector3D        v1;
  493.     TQ3Vector3D        v2;
  494.  
  495.     /* Search for the columns of the first two nonzero pivots */
  496.     for (j = 0; j < 3; j++) {
  497.         if (luDecomp->value[0][j] != 0.0) {
  498.             pivot0 = j;
  499.             break;
  500.         }
  501.     }
  502.  
  503.     if (pivot0 != -1) {
  504.         for (j = pivot0 + 1; j < 3; j++) {
  505.             if (luDecomp->value[1][j] != 0.0) {
  506.                 pivot1 = j;
  507.                 break;
  508.             }
  509.         }
  510.     }
  511.  
  512.     if (pivot0 == -1 || pivot1 == -1) {
  513.         /* Two nonzero pivots not found */
  514.         return kQ3Failure;
  515.     }
  516.  
  517.     /*
  518.      *  Points are represented as row vectors, and they are transformed as
  519.      *  follows:    p M = q
  520.      *  where M is a homogeneous matrix, and p and q are row vectors.
  521.      *  If M is the local-to-world matrix and A is the upper-left 3x3 submatrix
  522.      *  of M, then q is in the row space of A (in the absence of translation
  523.      *  and assuming that M is affine).  Given a factorization of A as
  524.      *  PA = LU, a basis for the row space of A is the basis of U.  Since
  525.      *  the rank of A is 2, the first two rows are a basis for U.  The space
  526.      *  spanned by these two rows is the plane in WC onto which all 3D geometry 
  527.      *  in LC is transformed.  The cross product of the two rows is the normal
  528.      *  of the plane.  Also note: the nullspace of A is orthogonal to the row 
  529.      *  space: this one dimensional space is spanned by the normal of the plane
  530.      *  in WC.
  531.      */
  532.     Q3Vector3D_Set(
  533.         &v0, 
  534.         luDecomp->value[0][0], 
  535.         luDecomp->value[0][1], 
  536.         luDecomp->value[0][2]);
  537.     Q3Vector3D_Set(&v1, 0.0, luDecomp->value[1][1], luDecomp->value[1][2]);
  538.     Q3Vector3D_Cross(&v0, &v1, &v2);
  539.     Q3Vector3D_Normalize(&v2, &flatWorld->normal);
  540.     
  541.     /*
  542.      *  A point on the plane in world coordinates (WC) is given by transforming
  543.      *  any point in local coordinates (LC) to WC.  If (0, 0, 0, 1) is the
  544.      *  LC point, the WC point is simply the last row of matrix M.  Substitute
  545.      *  this into the plane equation to obtain the constant D = - Ax - By - Cz.
  546.      */
  547.     flatWorld->constant = - matrix->value[3][0] * flatWorld->normal.x
  548.                           - matrix->value[3][1] * flatWorld->normal.y
  549.                           - matrix->value[3][2] * flatWorld->normal.z;
  550.     
  551.     /*
  552.      *  Each point in the column space of A has a one-to-one mapping to a point
  553.      *  in the row space of A.  In WC, the eye is on one side of the plane or
  554.      *  the other.  For the purposes of culling, it makes sense to put the eye
  555.      *  in LC on the same side of the plane.  So the eye vector in LC should be
  556.      *  orthogonal to the column space of A.  Otherwise, it would be in the 
  557.      *    column space, which doesn't make sense since the eye is almost never 
  558.      *  in the row space in WC.  A basis for the column space consists of the 
  559.      *  columns in A corresponding to those in U with nonzero pivots.
  560.      */
  561.     Q3Vector3D_Set(
  562.         &v0, 
  563.         matrix->value[0][pivot0],
  564.         matrix->value[1][pivot0], 
  565.         matrix->value[2][pivot0]);
  566.     Q3Vector3D_Set(
  567.         &v1, 
  568.         matrix->value[0][pivot1], 
  569.         matrix->value[1][pivot1], 
  570.         matrix->value[2][pivot1]);
  571.     Q3Vector3D_Cross(&v0, &v1, &v2);
  572.     Q3Vector3D_Normalize(&v2, &v2);
  573.     SRVector3D_To4D(&v2, eyeVectorLC);
  574.     
  575.     /*
  576.      *  We are almost done.  However, the sense of the eye vector in LC and
  577.      *  the normal of the plane equation in WC still need to be established.
  578.      *  The latter needs to point in the same hemisphere as the eye.
  579.      */
  580.     if (parallelProjection) {
  581.     
  582.         /* Parallel projection */
  583.         
  584.         SRVector4D_To3D(eyeVectorWC, &v2);
  585.         if (Q3Vector3D_Dot(&flatWorld->normal, &v2) < 0.0) {
  586.             /*
  587.              *  WC eye vector and plane normal point in opposite directions.
  588.              *  Need to reverse the sense of the plane equation.
  589.              */
  590.             Q3Vector3D_Negate(&flatWorld->normal, &flatWorld->normal);
  591.             flatWorld->constant = -flatWorld->constant;
  592.         }
  593.     } else {
  594.         TQ3Point3D    p0;
  595.         TQ3Point3D    p1;
  596.         
  597.         /*  Perspective projection */
  598.         
  599.         /*
  600.          *  A WC eye vector can be obtained by subtracting the eye point
  601.          *  by a point on the plane.  The LC point (0, 0, 0, 1) transformed
  602.          *  to WC is the last row of M.
  603.          */
  604.         Q3RationalPoint4D_To3D(eyePointWC, &p0);
  605.         Q3Point3D_Set(
  606.             &p1, 
  607.             matrix->value[3][0], 
  608.             matrix->value[3][1], 
  609.             matrix->value[3][2]);
  610.         Q3Point3D_Subtract(&p0, &p1, &v2);
  611.         if (Q3Vector3D_Dot(&flatWorld->normal, &v2) < 0.0) {
  612.             /*
  613.              *  WC eye vector and plane normal point in opposite directions.
  614.              *  Need to reverse the sense of the plane equation.
  615.              */
  616.             Q3Vector3D_Negate(&flatWorld->normal, &flatWorld->normal);
  617.             flatWorld->constant = -flatWorld->constant;
  618.         }
  619.     }
  620.     
  621.     /*
  622.      * The hard part to comprehend is the method for establishing the sense 
  623.      * of the eye vector in LC.  We had computed the eye vector by taking
  624.      * the cross product of v0 and v1, the two basis vectors for the column 
  625.      * space.  This gave a right-handed normal vector for the two basis 
  626.      * vectors.  If v0 and v1 are transformed by A, the two new vectors
  627.      * (v0 A) and (v1 A) lie in the row space, which is a plane parallel
  628.      * to the plane to which all geometry is transformed in WC.  The
  629.      * cross product of (v0 A) and (v1 A) either has the same or opposite
  630.      * sense as the normal to the plane, which points in the same hemisphere
  631.      * as the eye in WC.  If it is opposite, then the eye vector in LC has
  632.      * the opposite sense to the existing orientation.
  633.      */
  634.     Q3Vector3D_Transform(&v0, matrix, &v0);
  635.     Q3Vector3D_Transform(&v1, matrix, &v1);
  636.     Q3Vector3D_Cross(&v0, &v1, &v2);
  637.     if (Q3Vector3D_Dot(&flatWorld->normal, &v2) < 0.0) {
  638.         /*
  639.          * The LC eye vector transformed into WC and plane normal point in 
  640.          * opposite directions.  Need to reverse the sense of the LC eye
  641.          * vector.
  642.          */
  643.         SRVector4D_Negate(eyeVectorLC, eyeVectorLC);
  644.     }
  645.     
  646.     return kQ3Success;
  647. }
  648.  
  649.  
  650. /*===========================================================================*\
  651.  *
  652.  *    Routine:    SRVector4D_Set()
  653.  *
  654.  *    Comments:    Set a 4D vector
  655.  *
  656. \*===========================================================================*/
  657.  
  658. TSRVector4D *SRVector4D_Set(
  659.     TSRVector4D        *vector4D,
  660.     float            x, 
  661.     float            y,
  662.     float            z, 
  663.     float            w)
  664. {
  665.     assert(vector4D != NULL);
  666.     
  667.     vector4D->x = x;
  668.     vector4D->y = y;
  669.     vector4D->z = z;
  670.     vector4D->w = w;
  671.     
  672.     return (vector4D);
  673. }
  674.  
  675.  
  676. /*===========================================================================*\
  677.  *
  678.  *    Routine:    SRVector3D_To4D()
  679.  *
  680.  *    Comments:    Promote a 3D vector to 4D
  681.  *
  682. \*===========================================================================*/
  683.  
  684. TSRVector4D *SRVector3D_To4D(
  685.     const TQ3Vector3D    *vector3D,
  686.     TSRVector4D            *result)
  687. {
  688.     assert(vector3D != NULL);
  689.     assert(result != NULL);
  690.     
  691.     result->x = vector3D->x;
  692.     result->y = vector3D->y;
  693.     result->z = vector3D->z;
  694.     result->w = 1.0;
  695.     
  696.     return (result);
  697. }
  698.  
  699.  
  700. /*===========================================================================*\
  701.  *
  702.  *    Routine:    SRVector4D_To3D()
  703.  *
  704.  *    Comments:    Homogenize a 4D vector down to 3D
  705.  *
  706. \*===========================================================================*/
  707.  
  708. TQ3Vector3D *SRVector4D_To3D(
  709.     const TSRVector4D    *vector4D,
  710.     TQ3Vector3D            *result)
  711. {
  712.     float oneOverW;
  713.     
  714.     assert(vector4D != NULL);
  715.     assert(result != NULL);
  716.     
  717.     oneOverW = 1.0 / vector4D->w;
  718.     
  719.     result->x = vector4D->x * oneOverW;
  720.     result->y = vector4D->y * oneOverW;
  721.     result->z = vector4D->z * oneOverW;
  722.  
  723.     return (result);
  724. }
  725.  
  726.  
  727. /*===========================================================================*\
  728.  *
  729.  *    Routine:    SRVector4D_Negate()
  730.  *
  731.  *    Comments:    Negate a 4D vector
  732.  *
  733. \*===========================================================================*/
  734.  
  735. TSRVector4D *SRVector4D_Negate(
  736.     const TSRVector4D    *vector4D,
  737.     TSRVector4D            *result)
  738. {
  739.     assert(vector4D != NULL);
  740.     assert(result != NULL);
  741.  
  742.     result->x = -vector4D->x;
  743.     result->y = -vector4D->y;
  744.     result->z = -vector4D->z;
  745.     result->w = -vector4D->w;
  746.     
  747.     return (result);
  748. }
  749.  
  750.  
  751. /*===========================================================================*\
  752.  *
  753.  *    Routine:    SRVector4D_Normalize()
  754.  *
  755.  *    Comments:    Normalize a 4D vector
  756.  *
  757. \*===========================================================================*/
  758.  
  759. TSRVector4D *SRVector4D_Normalize(
  760.     const TSRVector4D    *vector4D,
  761.     TSRVector4D            *result)
  762. {
  763.     float    length, oneOverLength;
  764.     
  765.     assert(vector4D != NULL);
  766.     assert(result != NULL);
  767.  
  768.     length = sqrt(vector4D->x * vector4D->x + 
  769.                       vector4D->y * vector4D->y +
  770.                            vector4D->z * vector4D->z +
  771.                                vector4D->w * vector4D->w);
  772.     
  773.     /*
  774.      *  Check for zero-length vector
  775.      */
  776.     if (length < kQ3RealZero) {
  777.         *result = *vector4D;
  778.         return (NULL);
  779.     } else {
  780.         oneOverLength = 1.0 / length;
  781.     
  782.         result->x = vector4D->x * oneOverLength;
  783.         result->y = vector4D->y * oneOverLength;
  784.         result->z = vector4D->z * oneOverLength;
  785.         result->w = vector4D->w * oneOverLength;
  786.     }
  787.     
  788.     return (result);
  789. }
  790.  
  791.  
  792. /*===========================================================================*\
  793.  *
  794.  *    Routine:    SRVector4D_Transform()
  795.  *
  796.  *    Comments:    Transform a 4D vector by a matrix.
  797.  *
  798. \*===========================================================================*/
  799.     
  800. TSRVector4D *SRVector4D_Transform(
  801.     const TSRVector4D    *vector4D,
  802.     const TQ3Matrix4x4    *matrix4x4,
  803.     TSRVector4D            *result)
  804. {
  805.     TSRVector4D            s;
  806.     const TSRVector4D    *sPtr;
  807.     
  808.     assert(vector4D != NULL);
  809.     assert(matrix4x4 != NULL);
  810.     assert(result != NULL);
  811.  
  812.     /*
  813.      *  Check if caller is bashing the vector
  814.      */
  815.     if (vector4D == result) {
  816.         s = *vector4D;
  817.         sPtr = &s;
  818.     } else {
  819.         sPtr = vector4D;
  820.     }
  821.         
  822.     result->x = sPtr->x * matrix4x4->value[0][0] + 
  823.                 sPtr->y * matrix4x4->value[1][0] +
  824.                 sPtr->z * matrix4x4->value[2][0] +
  825.                 sPtr->w * matrix4x4->value[3][0];
  826.                 
  827.     result->y = sPtr->x * matrix4x4->value[0][1] +
  828.                 sPtr->y * matrix4x4->value[1][1] +
  829.                 sPtr->z * matrix4x4->value[2][1] +
  830.                 sPtr->w * matrix4x4->value[3][1];
  831.                 
  832.     result->z = sPtr->x * matrix4x4->value[0][2] + 
  833.                 sPtr->y * matrix4x4->value[1][2] +
  834.                 sPtr->z * matrix4x4->value[2][2] +
  835.                 sPtr->w * matrix4x4->value[3][2];
  836.                 
  837.     result->w = sPtr->x * matrix4x4->value[0][3] + 
  838.                 sPtr->y * matrix4x4->value[1][3] +
  839.                 sPtr->z * matrix4x4->value[2][3] +
  840.                 sPtr->w * matrix4x4->value[3][3];
  841.                 
  842.     return (result);
  843. }
  844.